Simulations results: different network topologies
Natasha Puzovic (puzovic@evolbio.mpg.de)
Nov 2, 2021
1 Read data
Read the genotype, phenotype and fitness information of evolved populations of 3000 network topology samples.
jointResultsFolder = "20211201_40node_0.05dens_joint_topos"
pathToJointResultsFolder = paste0("../results/", jointResultsFolder)
# read results and split into sel & neutr
allNetsResults_joint <- read.table(paste0(pathToJointResultsFolder, "/allNetsResults_prepped_joint.txt"),
sep = "\t", header = TRUE)
# rename nets
allNetsResults_joint[allNetsResults_joint$topo == "BA", "net"] <-
allNetsResults_joint[allNetsResults_joint$topo == "BA", "net"] + 1000
allNetsResults_joint[allNetsResults_joint$topo == "WS", "net"] <-
allNetsResults_joint[allNetsResults_joint$topo == "WS", "net"] + 2000
allNetsResults_joint$topo <- factor(allNetsResults_joint$topo, levels = c("ER", "BA", "WS"))
selResults <- allNetsResults_joint[allNetsResults_joint$scen == "sel", ]
neutrResults <- allNetsResults_joint[allNetsResults_joint$scen == "neu", ]
# subset responded genes
respondedToSelCutoff <- 0.5
# add which genes responded to selection
selResults$responseToSel <- selResults$s_g_area_abs > respondedToSelCutoff
# subset the genes that responded to selection
respondedGenes <- selResults[selResults$responseToSel == TRUE, ]1.1 Load packages, functions & colors
# Package list
packages <- c("nlme",
"MuMIn", # for r.squared in lme models
"ggplot2",
"ggpubr", # for the pubclean theme
"ggridges",
"gridExtra",
"cowplot", # for arranging plots
"infotheo",
"lme4",
"car", # for vif measure
"RColorBrewer", # for color palettes
"latex2exp", # for latex notation in the plots
"reshape2",
"Hmisc", # for correlation matrix
"corrplot", # for plotting correlation matrix
"ade4", # for PCA
"factoextra", # for scree plot
"FSA", # for Dunn tests
"rstatix", # for partial eta for lms
"dplyr", # for summarizing dataframes
"formattable", # for nice tables
"htmltools", # for outputting the table
"webshot") # for outputting the table
# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
# Packages loading
invisible(lapply(packages, library, character.only = TRUE))##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggpubr':
##
## get_legend## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmList## Loading required package: carData## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, units## corrplot 0.90 loaded## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa## Registered S3 methods overwritten by 'FSA':
## method from
## confint.boot car
## hist.boot car## ## FSA v0.9.3. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.##
## Attaching package: 'FSA'## The following object is masked from 'package:car':
##
## bootCase##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filter##
## Attaching package: 'dplyr'## The following objects are masked from 'package:Hmisc':
##
## src, summarize## The following object is masked from 'package:car':
##
## recode## The following object is masked from 'package:gridExtra':
##
## combine## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# MI
calcInformation <- function (v1, v2, binNum) {
# discretize
d_v1 <- discretize(v1, nbins = binNum);
d_v2 <- discretize(v2, nbins = binNum)
# mutual information
I_v1v2 <- mutinformation(d_v1, d_v2);
return("MI" = I_v1v2)
}
# colors for plots
noiseColor = "#01665e"
genotypeColor = "#7b3294"
phenotypeColor = "#d95f0e"
fitnessColor = "#e78ac3"
neutralityColor = "darkgray"
MIColor = "#2c7fb8"
topoColors = c("ER" = "darkgray", "BA" = "#c51b7d", "WS" = "#4d9221")netAllResults <- allNetsResults_joint %>%
group_by(scen, net) %>%
summarize(ave_varP_1 = mean(varP_1),
ave_varP_10000 = mean(varP_10000),
ave_relDeltaVar_10000 = mean(relDeltaVar_10000),
ave_s_g_area_abs = mean(s_g_area_abs),
topo = first(topo))## `summarise()` has grouped output by 'scen'. You can override using the `.groups` argument.netSelResults <- selResults %>%
group_by(scen, net) %>%
summarize(ave_varP_1 = mean(varP_1),
ave_varP_10000 = mean(varP_10000),
ave_relDeltaVar_10000 = mean(relDeltaVar_10000),
ave_s_g_area_abs = mean(s_g_area_abs),
ave_responseToSel = length(which(responseToSel)),
topo = first(topo))## `summarise()` has grouped output by 'scen'. You can override using the `.groups` argument.evolMetricsColnames <- c("meanG_1", "meanG_10000",
"meanP_1", "meanP_10000",
"varP_1", "varP_10000",
"CVP_1", "CVP_10000",
"noiseP_1", "noiseP_10000",
"FanoP_1", "FanoP_10000",
"relDeltaVar_1", "relDeltaVar_10000",
"relDeltaCV_1", "relDeltaCV_10000",
"relDeltaNoise_1", "relDeltaNoise_10000",
"relDeltaFano_1", "relDeltaFano_10000",
"s_g_area", "s_g_area_abs",
"s_p_area_relDeltaVar", "s_p_area_relDeltaCV",
"s_p_area_relDeltaNoise", "s_p_area_relDeltaFano")
evolMetricsColnames_ofInterest <- c("varP_1", "relDeltaVar_10000", "s_g_area_abs")
geneSpecificNetMetricsColnames <- c("k_all_inclps", "k_in_inclps", "k_out_inclps",
"clo_all", "betw", "eigen_centr",
"str_all_inclps", "str_in_inclps", "str_out_inclps",
"hub_score", "auth_incwght", "auth_excwght",
"absstr_all_inclps", "absInStr", "absOutStr",
"flow", "load", "info", "stress",
"absInStrT_sqrt", "absOutStrT_sqrt",
"absInStrT_log1p", "absOutStrT_log1p",
"absInStrT_log10", "absOutStrT_log10")
geneSpecificNetMetricsColnames_ofInterest <- c("k_all_inclps", "k_in_inclps", "k_out_inclps",
"clo_all", "betw", "eigen_centr",
"str_all_inclps", "str_in_inclps", "str_out_inclps",
"hub_score", "auth_incwght", "auth_excwght",
"absstr_all_inclps", "absInStr", "absOutStr",
"flow", "load", "info", "stress")
globalNetMetricsColnames <- c("diam", "meandst", "assort",
"cntr_degr_all", "cntr_indegr", "cntr_outdegr", "cntr_clo_all", "cntr_betw",
"ave_k_all_inclps", "ave_k_in_inclps","ave_k_out_inclps")
singletsSelResults <- selResults[!duplicated(selResults$net), c("net", globalNetMetricsColnames)]
netSelResults <- merge(netSelResults, singletsSelResults, by = "net")
singletsAllResults <- allNetsResults_joint[!duplicated(allNetsResults_joint$net), c("net", globalNetMetricsColnames)]
netAllResults <- merge(netAllResults, singletsAllResults, by = "net")2 Topological effects on distributions of noise metrics
2.1 Preview
# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = varP_1, x = topo)) +
geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = noiseColor) +
facet_wrap(~scen) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Expr. variance, gen 1"))) 
# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = varP_10000, x = topo)) +
geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = noiseColor) +
facet_wrap(~scen) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
ylab(expression(bold("Expr. variance, gen 10k"))) +
xlab("Network topology")
# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = relDeltaVar_10000, x = topo)) +
geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = noiseColor) +
facet_wrap(~scen) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
ylab(expression(bold("Rel."~Delta~"expr. variance"))) +
xlab("Network topology")## Warning: Removed 145 rows containing non-finite values (stat_boxplot).
# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = s_g_area_abs, x = topo)) +
geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = genotypeColor) +
facet_wrap(~scen) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = genotypeColor, "neutr" = genotypeColor)) +
ylab(expression(paste(bold("Selective pressure "), "|", bold(p), "|"))) +
xlab("Network topology")
my_comparisons <- list( c("ER", "BA"), c("BA", "WS"), c("WS", "ER") )
# just selection
plot_varFirstgen_sel <- ggplot(selResults, aes(y = varP_1, x = topo)) +
geom_violin(aes(fill = topo), trim = TRUE) +
geom_boxplot(width = 0.05, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label=round(..y.., digits=2))) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_fill_manual(values = topoColors) +
labs(x = "Network topology",
y = "Expr. variance, gen 1",
fill = "Topology")
plot_varFirstgen_sel
ggsave(filename = sprintf("plot_varFirstgen_sel.png"),
plot = plot_varFirstgen_sel,
path = pathToPlotsFolder,
device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
dpi = 300, limitsize = TRUE)
# just selection
plot_varLastgen_sel <- ggplot(selResults, aes(y = varP_10000, x = topo)) +
geom_violin(aes(fill = topo), trim = TRUE) +
geom_boxplot(width = 0.05, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label=round(..y.., digits=2))) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_fill_manual(values = topoColors) +
labs(x = "Network topology",
y = "Expr. variance, gen 10k",
fill = "Topology")
plot_varLastgen_sel
ggsave(filename = sprintf("plot_varLastgen_sel.png"),
plot = plot_varLastgen_sel,
path = pathToPlotsFolder,
device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
dpi = 300, limitsize = TRUE)
# just selection
plot_relDeltaVar_sel <- ggplot(selResults, aes(y = relDeltaVar_10000, x = topo)) +
geom_violin(aes(fill = topo), trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_fill_manual(values = topoColors) +
labs(x = "Network topology",
y = expression(bold("Rel."~Delta~"expr. variance")),
fill = "Topology")
plot_relDeltaVar_sel## Warning: Removed 101 rows containing non-finite values (stat_ydensity).## Warning: Removed 101 rows containing non-finite values (stat_boxplot).## Warning: Removed 101 rows containing non-finite values (stat_summary).## Warning: Removed 101 rows containing non-finite values (stat_signif).
ggsave(filename = sprintf("plot_relDeltaVar_sel.png"),
plot = plot_relDeltaVar_sel,
path = pathToPlotsFolder,
device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
dpi = 300, limitsize = TRUE)## Warning: Removed 101 rows containing non-finite values (stat_ydensity).## Warning: Removed 101 rows containing non-finite values (stat_boxplot).## Warning: Removed 101 rows containing non-finite values (stat_summary).## Warning: Removed 101 rows containing non-finite values (stat_signif).# just selection
plot_selpress_sel <- ggplot(selResults, aes(y = s_g_area_abs, x = topo)) +
geom_violin(aes(fill = topo), trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = 3, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_fill_manual(values = topoColors) +
labs(x = "Network topology",
y = expression(paste(bold("Selective pressure "), "|", bold(p), "|")),
fill = "Topology")
plot_selpress_sel
ggsave(filename = sprintf("plot_selpress_sel.png"),
plot = plot_selpress_sel,
path = pathToPlotsFolder,
device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
dpi = 300, limitsize = TRUE)
plot_violin_metrics_per_topo <- plot_grid(plot_varFirstgen_sel,
plot_varLastgen_sel,
plot_relDeltaVar_sel,
plot_selpress_sel,
scale = c(0.95, 0.95, 0.95, 0.95),
labels = "AUTO",
label_size = 20,
label_fontface = "bold")## Warning: Removed 101 rows containing non-finite values (stat_ydensity).## Warning: Removed 101 rows containing non-finite values (stat_boxplot).## Warning: Removed 101 rows containing non-finite values (stat_summary).## Warning: Removed 101 rows containing non-finite values (stat_signif).ggsave(filename = sprintf("plot_violin_metrics_per_topo.png"),
plot = plot_violin_metrics_per_topo,
bg = "white",
path = pathToPlotsFolder,
device = "png", scale = 1.8, width = 17, height = 12, units = "cm",
dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_violin_metrics_per_topo.tiff"),
plot = plot_violin_metrics_per_topo,
bg = "white",
path = pathToPlotsFolder,
device = "tiff", scale = 1.8, width = 17, height = 12, units = "cm",
dpi = 300, limitsize = TRUE)2.2 Comparisons
summary(selResults$varP_1[selResults$topo == "BA"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0869 187.5738 323.5003 381.4193 505.1286 2166.1780summary(selResults$varP_1[selResults$topo == "ER"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0001 97.7791 198.1725 325.6116 391.2855 2397.0850summary(selResults$varP_1[selResults$topo == "WS"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0001 117.6926 231.2422 389.6922 501.9861 2326.2150summary(selResults$varP_10000[selResults$topo == "BA"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 8.587 12.434 14.055 16.863 572.432summary(selResults$varP_10000[selResults$topo == "ER"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 3.658 8.220 20.781 15.763 2425.164summary(selResults$varP_10000[selResults$topo == "WS"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 3.301 8.052 24.449 18.659 2130.139summary(selResults$relDeltaVar_10000[selResults$topo == "BA"])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1.00000 -0.94763 -0.93174 -0.92602 -0.90960 0.01953 9summary(selResults$relDeltaVar_10000[selResults$topo == "ER"])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1.000 -0.951 -0.926 -0.911 -0.894 1.000 37summary(selResults$relDeltaVar_10000[selResults$topo == "WS"])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1.0000 -0.9591 -0.9357 -0.9172 -0.9046 0.5278 55summary(selResults$s_g_area_abs[selResults$topo == "BA"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000025 0.8314451 0.8538640 0.7917163 0.8762081 0.9862562summary(selResults$s_g_area_abs[selResults$topo == "ER"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00001 0.82969 0.88379 0.78503 0.93096 0.98452summary(selResults$s_g_area_abs[selResults$topo == "WS"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000021 0.7876157 0.8919116 0.7460259 0.9367470 0.9858500kruskal.test(varP_1 ~ topo, data = selResults)##
## Kruskal-Wallis rank sum test
##
## data: varP_1 by topo
## Kruskal-Wallis chi-squared = 3118.6, df = 2, p-value < 2.2e-16kruskal.test(varP_10000 ~ topo, data = selResults)##
## Kruskal-Wallis rank sum test
##
## data: varP_10000 by topo
## Kruskal-Wallis chi-squared = 2620.6, df = 2, p-value < 2.2e-16kruskal.test(relDeltaVar_10000 ~ topo, data = selResults)##
## Kruskal-Wallis rank sum test
##
## data: relDeltaVar_10000 by topo
## Kruskal-Wallis chi-squared = 900.53, df = 2, p-value < 2.2e-16kruskal.test(s_g_area_abs ~ topo, data = selResults)##
## Kruskal-Wallis rank sum test
##
## data: s_g_area_abs by topo
## Kruskal-Wallis chi-squared = 3428, df = 2, p-value < 2.2e-16pairwise.wilcox.test(selResults$varP_1, selResults$topo, paired = FALSE, p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: selResults$varP_1 and selResults$topo
##
## ER BA
## BA <2e-16 -
## WS <2e-16 <2e-16
##
## P value adjustment method: BHpairwise.wilcox.test(selResults$varP_10000, selResults$topo, paired = FALSE, p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: selResults$varP_10000 and selResults$topo
##
## ER BA
## BA <2e-16 -
## WS 0.86 <2e-16
##
## P value adjustment method: BHpairwise.wilcox.test(selResults$relDeltaVar_10000, selResults$topo, paired = FALSE, p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: selResults$relDeltaVar_10000 and selResults$topo
##
## ER BA
## BA <2e-16 -
## WS <2e-16 <2e-16
##
## P value adjustment method: BHpairwise.wilcox.test(selResults$s_g_area_abs, selResults$topo, paired = FALSE, p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: selResults$s_g_area_abs and selResults$topo
##
## ER BA
## BA <2e-16 -
## WS 0.064 <2e-16
##
## P value adjustment method: BH# ER to BA
wilcox.test(selResults$varP_1[selResults$topo == "ER"], selResults$varP_1[selResults$topo == "BA"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_1[selResults$topo == "ER"] and selResults$varP_1[selResults$topo == "BA"]
## W = 559264846, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$varP_1[selResults$topo == "ER"], selResults$varP_1[selResults$topo == "BA"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_1[selResults$topo == "ER"] and selResults$varP_1[selResults$topo == "BA"]
## W = 559264846, p-value = 1
## alternative hypothesis: true location shift is greater than 0# BA to WS
wilcox.test(selResults$varP_1[selResults$topo == "BA"], selResults$varP_1[selResults$topo == "WS"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_1[selResults$topo == "BA"] and selResults$varP_1[selResults$topo == "WS"]
## W = 816760039, p-value = 1
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$varP_1[selResults$topo == "BA"], selResults$varP_1[selResults$topo == "WS"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_1[selResults$topo == "BA"] and selResults$varP_1[selResults$topo == "WS"]
## W = 816760039, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0# WS to ER
wilcox.test(selResults$varP_1[selResults$topo == "WS"], selResults$varP_1[selResults$topo == "ER"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_1[selResults$topo == "WS"] and selResults$varP_1[selResults$topo == "ER"]
## W = 745782968, p-value = 1
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$varP_1[selResults$topo == "WS"], selResults$varP_1[selResults$topo == "ER"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_1[selResults$topo == "WS"] and selResults$varP_1[selResults$topo == "ER"]
## W = 745782968, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0# ER to BA
wilcox.test(selResults$varP_10000[selResults$topo == "ER"], selResults$varP_10000[selResults$topo == "BA"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_10000[selResults$topo == "ER"] and selResults$varP_10000[selResults$topo == "BA"]
## W = 586187816, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$varP_10000[selResults$topo == "ER"], selResults$varP_10000[selResults$topo == "BA"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_10000[selResults$topo == "ER"] and selResults$varP_10000[selResults$topo == "BA"]
## W = 586187816, p-value = 1
## alternative hypothesis: true location shift is greater than 0# BA to WS
wilcox.test(selResults$varP_10000[selResults$topo == "BA"], selResults$varP_10000[selResults$topo == "WS"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_10000[selResults$topo == "BA"] and selResults$varP_10000[selResults$topo == "WS"]
## W = 843433923, p-value = 1
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$varP_10000[selResults$topo == "BA"], selResults$varP_10000[selResults$topo == "WS"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_10000[selResults$topo == "BA"] and selResults$varP_10000[selResults$topo == "WS"]
## W = 843433923, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0# WS to ER
wilcox.test(selResults$varP_10000[selResults$topo == "WS"], selResults$varP_10000[selResults$topo == "ER"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_10000[selResults$topo == "WS"] and selResults$varP_10000[selResults$topo == "ER"]
## W = 683677545, p-value = 0.4305
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$varP_10000[selResults$topo == "WS"], selResults$varP_10000[selResults$topo == "ER"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$varP_10000[selResults$topo == "WS"] and selResults$varP_10000[selResults$topo == "ER"]
## W = 683677545, p-value = 0.5695
## alternative hypothesis: true location shift is greater than 0# ER to BA
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "ER"], selResults$relDeltaVar_10000[selResults$topo == "BA"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$relDeltaVar_10000[selResults$topo == "ER"] and selResults$relDeltaVar_10000[selResults$topo == "BA"]
## W = 789419244, p-value = 1
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "ER"], selResults$relDeltaVar_10000[selResults$topo == "BA"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$relDeltaVar_10000[selResults$topo == "ER"] and selResults$relDeltaVar_10000[selResults$topo == "BA"]
## W = 789419244, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0# BA to WS
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "BA"], selResults$relDeltaVar_10000[selResults$topo == "WS"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$relDeltaVar_10000[selResults$topo == "BA"] and selResults$relDeltaVar_10000[selResults$topo == "WS"]
## W = 762636206, p-value = 1
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "BA"], selResults$relDeltaVar_10000[selResults$topo == "WS"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$relDeltaVar_10000[selResults$topo == "BA"] and selResults$relDeltaVar_10000[selResults$topo == "WS"]
## W = 762636206, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0# WS to ER
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "WS"], selResults$relDeltaVar_10000[selResults$topo == "ER"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$relDeltaVar_10000[selResults$topo == "WS"] and selResults$relDeltaVar_10000[selResults$topo == "ER"]
## W = 600249666, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "WS"], selResults$relDeltaVar_10000[selResults$topo == "ER"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$relDeltaVar_10000[selResults$topo == "WS"] and selResults$relDeltaVar_10000[selResults$topo == "ER"]
## W = 600249666, p-value = 1
## alternative hypothesis: true location shift is greater than 0# ER to BA
wilcox.test(selResults$s_g_area_abs[selResults$topo == "ER"], selResults$s_g_area_abs[selResults$topo == "BA"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$s_g_area_abs[selResults$topo == "ER"] and selResults$s_g_area_abs[selResults$topo == "BA"]
## W = 895602570, p-value = 1
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$s_g_area_abs[selResults$topo == "ER"], selResults$s_g_area_abs[selResults$topo == "BA"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$s_g_area_abs[selResults$topo == "ER"] and selResults$s_g_area_abs[selResults$topo == "BA"]
## W = 895602570, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0# BA to WS
wilcox.test(selResults$s_g_area_abs[selResults$topo == "BA"], selResults$s_g_area_abs[selResults$topo == "WS"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$s_g_area_abs[selResults$topo == "BA"] and selResults$s_g_area_abs[selResults$topo == "WS"]
## W = 578313394, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$s_g_area_abs[selResults$topo == "BA"], selResults$s_g_area_abs[selResults$topo == "WS"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$s_g_area_abs[selResults$topo == "BA"] and selResults$s_g_area_abs[selResults$topo == "WS"]
## W = 578313394, p-value = 1
## alternative hypothesis: true location shift is greater than 0# WS to ER
wilcox.test(selResults$s_g_area_abs[selResults$topo == "WS"], selResults$s_g_area_abs[selResults$topo == "ER"],
alternative = "less")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$s_g_area_abs[selResults$topo == "WS"] and selResults$s_g_area_abs[selResults$topo == "ER"]
## W = 689574788, p-value = 0.9682
## alternative hypothesis: true location shift is less than 0wilcox.test(selResults$s_g_area_abs[selResults$topo == "WS"], selResults$s_g_area_abs[selResults$topo == "ER"],
alternative = "greater")##
## Wilcoxon rank sum test with continuity correction
##
## data: selResults$s_g_area_abs[selResults$topo == "WS"] and selResults$s_g_area_abs[selResults$topo == "ER"]
## W = 689574788, p-value = 0.03178
## alternative hypothesis: true location shift is greater than 0dunnTest(varP_1 ~ topo, data = selResults, method = "bh")## Dunn (1964) Kruskal-Wallis multiple comparison## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BA - ER 55.45841 0.000000e+00 0.000000e+00
## 2 BA - WS 32.62836 1.624895e-233 2.437343e-233
## 3 ER - WS -22.33020 1.880903e-110 1.880903e-110dunnTest(varP_10000 ~ topo, data = selResults, method = "bh")## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BA - ER 45.327324 0.0000000 0.0000000
## 2 BA - WS 42.787665 0.0000000 0.0000000
## 3 ER - WS -2.338746 0.0193486 0.0193486dunnTest(relDeltaVar_10000 ~ topo, data = selResults, method = "bh")## Warning: Some rows deleted from 'x' and 'g' because missing data.## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BA - ER -17.67761 6.238093e-70 9.357140e-70
## 2 BA - WS 12.64267 1.228083e-36 1.228083e-36
## 3 ER - WS 29.85645 7.240590e-196 2.172177e-195dunnTest(s_g_area_abs ~ topo, data = selResults, method = "bh")## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BA - ER -50.9647916 0.0000000 0.0000000
## 2 BA - WS -49.9192228 0.0000000 0.0000000
## 3 ER - WS 0.8429799 0.3992397 0.3992397lmeModel <- lm(ave_varP_1 ~ topo,
data = netSelResults)
summary(lmeModel)##
## Call:
## lm(formula = ave_varP_1 ~ topo, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -295.62 -71.22 -10.35 58.69 514.02
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 324.224 3.190 101.62 <2e-16 ***
## topoBA 55.686 4.512 12.34 <2e-16 ***
## topoWS 63.437 4.512 14.06 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 100.9 on 2997 degrees of freedom
## Multiple R-squared: 0.0728, Adjusted R-squared: 0.07218
## F-statistic: 117.6 on 2 and 2997 DF, p-value: < 2.2e-16r.squaredGLMM(lmeModel)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.07275108 0.07275108lmeModel <- lm(ave_varP_10000 ~ topo,
data = netSelResults)
summary(lmeModel)##
## Call:
## lm(formula = ave_varP_10000 ~ topo, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.04 -8.98 -3.34 0.36 762.49
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.5092 0.9596 21.373 < 2e-16 ***
## topoBA -6.5096 1.3571 -4.797 1.69e-06 ***
## topoWS 3.5666 1.3571 2.628 0.00863 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.35 on 2997 degrees of freedom
## Multiple R-squared: 0.01857, Adjusted R-squared: 0.01791
## F-statistic: 28.35 on 2 and 2997 DF, p-value: 6.361e-13r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.01855435 0.01855435lmeModel <- lm(ave_relDeltaVar_10000 ~ topo,
data = netSelResults[!is.na(netSelResults$ave_relDeltaVar_10000), ])
summary(lmeModel)##
## Call:
## lm(formula = ave_relDeltaVar_10000 ~ topo, data = netSelResults[!is.na(netSelResults$ave_relDeltaVar_10000),
## ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03773 -0.01319 -0.00588 0.00243 0.60125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.911215 0.001181 -771.311 < 2e-16 ***
## topoBA -0.014871 0.001660 -8.958 < 2e-16 ***
## topoWS -0.006834 0.001680 -4.069 4.85e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03672 on 2900 degrees of freedom
## Multiple R-squared: 0.027, Adjusted R-squared: 0.02633
## F-statistic: 40.23 on 2 and 2900 DF, p-value: < 2.2e-16r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.02698037 0.02698037lmeModel <- lm(ave_s_g_area_abs ~ topo,
data = netSelResults)
summary(lmeModel)##
## Call:
## lm(formula = ave_s_g_area_abs ~ topo, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63787 -0.02335 0.01166 0.03866 0.12709
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.785488 0.002050 383.210 <2e-16 ***
## topoBA 0.004880 0.002899 1.683 0.0924 .
## topoWS -0.038750 0.002899 -13.367 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06482 on 2997 degrees of freedom
## Multiple R-squared: 0.08322, Adjusted R-squared: 0.0826
## F-statistic: 136 on 2 and 2997 DF, p-value: < 2.2e-16r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.08316459 0.083164593 Topological effects on local network metrics
3.1 Logistic regression
logRegModel <- glmer(responseToSel ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo + (1|net),
data = selResults,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
summary(logRegModel)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: responseToSel ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo + (1 |
## net)
## Data: selResults
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 69700.7 69797.0 -34840.3 69680.7 113264
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.8288 0.1336 0.2236 0.3491 5.5940
##
## Random effects:
## Groups Name Variance Std.Dev.
## net (Intercept) 0.3268 0.5717
## Number of obs: 113274, groups: net, 3000
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.1053574 0.0645610 79.078 < 2e-16 ***
## absInStrT_sqrt -1.9270045 0.0284020 -67.848 < 2e-16 ***
## absOutStrT_sqrt -0.0829021 0.0226893 -3.654 0.000258 ***
## topoBA 0.9209796 0.1075134 8.566 < 2e-16 ***
## topoWS -0.2684890 0.0945709 -2.839 0.004525 **
## absInStrT_sqrt:topoBA 0.0120205 0.0516509 0.233 0.815975
## absInStrT_sqrt:topoWS 0.0006353 0.0401797 0.016 0.987384
## absOutStrT_sqrt:topoBA -0.2947402 0.0252770 -11.660 < 2e-16 ***
## absOutStrT_sqrt:topoWS -0.0728901 0.0333220 -2.187 0.028710 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W aOST_:B
## absInStrT_s -0.819
## absOtStrT_s -0.464 0.064
## topoBA -0.584 0.477 0.280
## topoWS -0.677 0.554 0.317 0.403
## absInST_:BA 0.438 -0.540 -0.036 -0.899 -0.302
## absInST_:WS 0.573 -0.702 -0.045 -0.341 -0.818 0.384
## absOtST_:BA 0.412 -0.053 -0.898 -0.384 -0.284 0.132 0.040
## absOtST_:WS 0.316 -0.044 -0.681 -0.190 -0.560 0.024 0.154 0.6113.2 Local network metrics
lmeModel <- lme(varP_1 ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo,
data = selResults,
weights = varExp(form = ~absInStrT_sqrt),
random = ~1|net,
method = "ML")
summary(lmeModel)## Linear mixed-effects model fit by maximum likelihood
## Data: selResults
## AIC BIC logLik
## 1535840 1535956 -767908.1
##
## Random effects:
## Formula: ~1 | net
## (Intercept) Residual
## StdDev: 17.62603 43.53776
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~absInStrT_sqrt
## Parameter estimates:
## expon
## 1.160825
## Fixed effects: varP_1 ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo
## Value Std.Error DF t-value p-value
## (Intercept) 98.39128 1.4085297 110268 69.85389 0
## absInStrT_sqrt 160.81561 1.1328027 110268 141.96259 0
## absOutStrT_sqrt -6.24249 0.7220382 110268 -8.64565 0
## topoBA -15.23558 2.8368173 2997 -5.37066 0
## topoWS -31.08754 2.9511842 2997 -10.53392 0
## absInStrT_sqrt:topoBA 13.00962 2.1763677 110268 5.97768 0
## absInStrT_sqrt:topoWS 41.02790 1.8236006 110268 22.49829 0
## absOutStrT_sqrt:topoBA 3.97398 0.8717507 110268 4.55862 0
## absOutStrT_sqrt:topoWS 8.91289 1.3415885 110268 6.64353 0
## Correlation:
## (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W
## absInStrT_sqrt -0.430
## absOutStrT_sqrt -0.823 0.267
## topoBA -0.497 0.214 0.409
## topoWS -0.477 0.205 0.393 0.237
## absInStrT_sqrt:topoBA 0.224 -0.521 -0.139 -0.744 -0.107
## absInStrT_sqrt:topoWS 0.267 -0.621 -0.166 -0.133 -0.579 0.323
## absOutStrT_sqrt:topoBA 0.682 -0.221 -0.828 -0.741 -0.325 0.453 0.138
## absOutStrT_sqrt:topoWS 0.443 -0.144 -0.538 -0.220 -0.884 0.075 0.402
## aOST_:B
## absInStrT_sqrt
## absOutStrT_sqrt
## topoBA
## topoWS
## absInStrT_sqrt:topoBA
## absInStrT_sqrt:topoWS
## absOutStrT_sqrt:topoBA
## absOutStrT_sqrt:topoWS 0.446
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.23168986 -0.60002511 -0.07796252 0.46211907 9.18505792
##
## Number of Observations: 113274
## Number of Groups: 3000r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.8665274 0.8853229lmeModel <- lme(relDeltaVar_10000 ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo,
data = respondedGenes,
weights = varExp(form = ~absInStrT_sqrt),
random = ~1|net,
method = "ML")
summary(lmeModel)## Linear mixed-effects model fit by maximum likelihood
## Data: respondedGenes
## AIC BIC logLik
## -315994.1 -315880.1 158009
##
## Random effects:
## Formula: ~1 | net
## (Intercept) Residual
## StdDev: 0.01993878 0.03603974
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~absInStrT_sqrt
## Parameter estimates:
## expon
## 0.2132598
## Fixed effects: relDeltaVar_10000 ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo
## Value Std.Error DF t-value p-value
## (Intercept) -0.8892058 0.0009009564 95608 -986.9577 0
## absInStrT_sqrt 0.0078846 0.0003751123 95608 21.0193 0
## absOutStrT_sqrt -0.0236606 0.0003348018 95608 -70.6705 0
## topoBA 0.0063155 0.0014154920 2997 4.4617 0
## topoWS -0.0157462 0.0015008533 2997 -10.4915 0
## absInStrT_sqrt:topoBA -0.0282349 0.0006710705 95608 -42.0743 0
## absInStrT_sqrt:topoWS 0.0051777 0.0006073184 95608 8.5256 0
## absOutStrT_sqrt:topoBA 0.0084510 0.0003862807 95608 21.8778 0
## absOutStrT_sqrt:topoWS 0.0033987 0.0005645936 95608 6.0198 0
## Correlation:
## (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W
## absInStrT_sqrt -0.467
## absOutStrT_sqrt -0.561 0.228
## topoBA -0.636 0.297 0.357
## topoWS -0.600 0.280 0.337 0.382
## absInStrT_sqrt:topoBA 0.261 -0.559 -0.127 -0.660 -0.157
## absInStrT_sqrt:topoWS 0.288 -0.618 -0.141 -0.184 -0.568 0.345
## absOutStrT_sqrt:topoBA 0.486 -0.198 -0.867 -0.491 -0.292 0.308 0.122
## absOutStrT_sqrt:topoWS 0.333 -0.135 -0.593 -0.212 -0.661 0.076 0.306
## aOST_:B
## absInStrT_sqrt
## absOutStrT_sqrt
## topoBA
## topoWS
## absInStrT_sqrt:topoBA
## absInStrT_sqrt:topoWS
## absOutStrT_sqrt:topoBA
## absOutStrT_sqrt:topoWS 0.514
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -8.0758560 -0.4668731 -0.1549897 0.2764916 30.5878268
##
## Number of Observations: 98614
## Number of Groups: 3000r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.1832294 0.3746393lmeModel <- lme(s_g_area_abs ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo,
data = respondedGenes,
weights = varExp(form = ~absInStrT_sqrt),
random = ~1|net,
method = "ML")
summary(lmeModel)## Linear mixed-effects model fit by maximum likelihood
## Data: respondedGenes
## AIC BIC logLik
## -298672.5 -298558.5 149348.2
##
## Random effects:
## Formula: ~1 | net
## (Intercept) Residual
## StdDev: 0.01277971 0.02940703
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~absInStrT_sqrt
## Parameter estimates:
## expon
## 0.450376
## Fixed effects: s_g_area_abs ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo
## Value Std.Error DF t-value p-value
## (Intercept) 0.8823365 0.0007457204 95608 1183.2001 0.0000
## absInStrT_sqrt -0.0376771 0.0004010174 95608 -93.9538 0.0000
## absOutStrT_sqrt 0.0347294 0.0003346249 95608 103.7860 0.0000
## topoBA 0.0018673 0.0012663086 2997 1.4746 0.1404
## topoWS 0.0222238 0.0013542170 2997 16.4108 0.0000
## absInStrT_sqrt:topoBA 0.0143136 0.0007355510 95608 19.4598 0.0000
## absInStrT_sqrt:topoWS -0.0055371 0.0006547984 95608 -8.4562 0.0000
## absOutStrT_sqrt:topoBA -0.0151242 0.0003887359 95608 -38.9061 0.0000
## absOutStrT_sqrt:topoWS -0.0075333 0.0005780342 95608 -13.0327 0.0000
## Correlation:
## (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W
## absInStrT_sqrt -0.513
## absOutStrT_sqrt -0.696 0.261
## topoBA -0.589 0.302 0.410
## topoWS -0.551 0.283 0.383 0.324
## absInStrT_sqrt:topoBA 0.280 -0.545 -0.142 -0.747 -0.154
## absInStrT_sqrt:topoWS 0.314 -0.612 -0.160 -0.185 -0.619 0.334
## absOutStrT_sqrt:topoBA 0.599 -0.224 -0.861 -0.609 -0.330 0.368 0.137
## absOutStrT_sqrt:topoWS 0.403 -0.151 -0.579 -0.237 -0.777 0.082 0.347
## aOST_:B
## absInStrT_sqrt
## absOutStrT_sqrt
## topoBA
## topoWS
## absInStrT_sqrt:topoBA
## absInStrT_sqrt:topoWS
## absOutStrT_sqrt:topoBA
## absOutStrT_sqrt:topoWS 0.498
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -10.2408639 -0.3200772 0.1428952 0.5748408 3.5412652
##
## Number of Observations: 98614
## Number of Groups: 3000r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.5930888 0.6577299lmeModel <- lme(s_g_area_abs ~ absInStrT_sqrt + absOutStrT_sqrt + topo,
data = respondedGenes,
weights = varExp(form = ~absInStrT_sqrt),
random = ~1|net,
method = "ML")
summary(lmeModel)## Linear mixed-effects model fit by maximum likelihood
## Data: respondedGenes
## AIC BIC logLik
## -295317.6 -295241.7 147666.8
##
## Random effects:
## Formula: ~1 | net
## (Intercept) Residual
## StdDev: 0.01241693 0.03028914
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~absInStrT_sqrt
## Parameter estimates:
## expon
## 0.4413942
## Fixed effects: s_g_area_abs ~ absInStrT_sqrt + absOutStrT_sqrt + topo
## Value Std.Error DF t-value p-value
## (Intercept) 0.8997841 0.0005686180 95612 1582.4053 0.0000
## absInStrT_sqrt -0.0372239 0.0002857357 95612 -130.2740 0.0000
## absOutStrT_sqrt 0.0215142 0.0001524232 95612 141.1476 0.0000
## topoBA -0.0001470 0.0006926962 2997 -0.2123 0.8319
## topoWS 0.0080884 0.0006752761 2997 11.9779 0.0000
## Correlation:
## (Intr) abIST_ abOST_ topoBA
## absInStrT_sqrt -0.476
## absOutStrT_sqrt -0.491 0.390
## topoBA -0.470 -0.208 -0.014
## topoWS -0.484 -0.126 -0.090 0.482
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -10.0788670 -0.3343392 0.1765577 0.5932647 3.5587811
##
## Number of Observations: 98614
## Number of Groups: 3000r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.580069 0.64048744 Correlograms
# subset
ave_evolMetricsColnames_ofInterest <- c("ave_varP_1", "ave_relDeltaVar_10000", "ave_s_g_area_abs")
# There are also "num_generations", "num_nodes", "dens", "pop_size" columns, but they are identical for all topos.
# global net metrics
globalMetrics_forCorrs_sel <- netSelResults[, c(ave_evolMetricsColnames_ofInterest, globalNetMetricsColnames)]
globalMetrics_forCorrs_neu <- netSelResults[, c(ave_evolMetricsColnames_ofInterest, globalNetMetricsColnames)]
# gene-specific metrics
geneMetrics_forCorrs_sel <- selResults[, c(evolMetricsColnames_ofInterest, geneSpecificNetMetricsColnames_ofInterest)]
geneMetrics_forCorrs_neu <- neutrResults[, c(evolMetricsColnames_ofInterest, geneSpecificNetMetricsColnames_ofInterest)]4.1 Gene-specific metrics
# selection
geneSpecificCorrs_sel <- rcorr(as.matrix(geneMetrics_forCorrs_sel), type = "spearman")
# neutrality
geneSpecificCorrs_neu <- rcorr(as.matrix(geneMetrics_forCorrs_neu), type = "spearman")png(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_sel.png"),
width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.7,
p.mat = geneSpecificCorrs_sel$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 2tiff(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_sel.tiff"),
width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.7,
p.mat = geneSpecificCorrs_sel$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 2png(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_neu.png"),
width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.7,
p.mat = geneSpecificCorrs_neu$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 2tiff(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_neu.tiff"),
width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.7,
p.mat = geneSpecificCorrs_neu$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 24.2 Global network metrics
# selection
globalCorrs_sel <- rcorr(as.matrix(globalMetrics_forCorrs_sel), type = "spearman")
# neutrality
globalCorrs_neu <- rcorr(as.matrix(globalMetrics_forCorrs_neu), type = "spearman")png(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_sel.png"),
width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.8,
p.mat = globalCorrs_sel$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 2tiff(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_sel.tiff"),
width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.8,
p.mat = globalCorrs_sel$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 2png(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_neu.png"),
width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.8,
p.mat = globalCorrs_neu$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 2tiff(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_neu.tiff"),
width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
number.cex = 0.8,
p.mat = globalCorrs_neu$P, sig.level = 0.05, insig = "blank",
tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"), outline = TRUE, mar=c(0,0,1,0))
dev.off()## png
## 25 PCA of global network metrics
globalMetrics_forPCA <- netAllResults[, globalNetMetricsColnames]
colnames(globalMetrics_forPCA) <- c("diameter", "mean path distance", "degree assortativity", "degree centralization",
"indegree centralization", "outdegree centralization", "closeness centralization",
"betweenness centralization", "average degree", "average indegree", "average outdegree")
dudi <- dudi.pca(globalMetrics_forPCA, center = TRUE, scale = TRUE, nf = 10, scannf = FALSE)
plot_biplot <- fviz_pca_biplot(dudi,
geom.ind = "point",
col.ind = netAllResults$topo,
col.var = "black",
repel = TRUE,
addEllipses = TRUE,
legend.title = "Topology") +
scale_colour_manual(values = topoColors)
plot_biplot
plot_scree <- fviz_eig(dudi, addlabels = TRUE)
plot_corCircle <- fviz_pca_var(dudi, col.var = "contrib", labelsize = 4, repel = TRUE) +
scale_color_gradient2(low = "white", mid = "blue", high = "black")
plot_PCA <- plot_grid(plot_scree, plot_corCircle,
scale = c(0.95, 0.95),
labels = "AUTO",
label_size = 20,
label_fontface = "bold",
ncol = 2)
ggsave(filename = sprintf("plot_PCA.png"),
plot = plot_PCA,
path = pathToPlotsFolder, bg = 'white',
device = "png", scale = 2, width = 18, height = 9, units = "cm",
dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_PCA.tiff"),
plot = plot_PCA,
path = pathToPlotsFolder, bg = 'white',
device = "tiff", scale = 2, width = 18, height = 9, units = "cm",
dpi = 300, limitsize = TRUE)
netAllResults$PC1<- -dudi$li$Axis1
netAllResults$PC2<- -dudi$li$Axis2
ave_responseToSel <- netSelResults$ave_responseToSel
netSelResults <- netAllResults[netAllResults$scen == "sel", ]
netSelResults$ave_responseToSel <- ave_responseToSel
netNeuResults <- netAllResults[netAllResults$scen == "neu", ]# just selection
plot_PC1_topos <- ggplot(netSelResults, aes(y = PC1, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("PC1 (diameter + centralization)")))
plot_PC1_topos
plot_PC2_topos <- ggplot(netSelResults, aes(y = PC2, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("PC2 (average degree)")))
plot_PC2_topos
plot_netMetric <- ggplot(netSelResults, aes(y = diam, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Diameter")))
plot_netMetric
plot_netMetric <- ggplot(netSelResults, aes(y = cntr_indegr, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Indegree centralization")))
plot_netMetric
plot_netMetric <- ggplot(netSelResults, aes(y = cntr_outdegr, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Outdegree centralization")))
plot_netMetric
plot_netMetric <- ggplot(netSelResults, aes(y = ave_k_all_inclps, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Average degree")))
plot_netMetric
plot_netMetric <- ggplot(netSelResults, aes(y = ave_k_in_inclps, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Average indegree")))
plot_netMetric
plot_netMetric <- ggplot(netSelResults, aes(y = ave_k_out_inclps, x = topo)) +
geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
labs(x = "Network topology",
y = expression(bold("Average outdegree")))
plot_netMetric
6 LMs with PCs
6.1 Number of responded genes
lmModel <- lm(ave_responseToSel ~ PC1 + PC2,
data = netSelResults)
summary(lmModel)##
## Call:
## lm(formula = ave_responseToSel ~ PC1 + PC2, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.9806 -1.2277 0.5615 2.0439 7.7641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.87133 0.06275 523.82 <2e-16 ***
## PC1 -0.71303 0.02521 -28.29 <2e-16 ***
## PC2 -0.76930 0.03506 -21.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.437 on 2997 degrees of freedom
## Multiple R-squared: 0.2995, Adjusted R-squared: 0.2991
## F-statistic: 640.8 on 2 and 2997 DF, p-value: < 2.2e-16plot_PC1 <- ggplot(netSelResults, aes(y = ave_responseToSel, x = PC1)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
geom_abline(slope = lmModel$coefficients[2],
intercept = lmModel$coefficients[1],
color = "black", linetype = "dashed", size = 1) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_colour_manual(values = topoColors) +
labs(x = expression(bold("PC1 (diameter + centralization)")),
y = expression(paste(bold("Average # responded genes"))))
plot_PC1
plot_PC2 <- ggplot(netSelResults, aes(y = ave_responseToSel, x = PC2)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
geom_abline(slope = lmModel$coefficients[3],
intercept = lmModel$coefficients[1],
color = "black", linetype = "dashed", size = 1) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_colour_manual(values = topoColors) +
labs(x = expression(bold("PC2 (average degree)")),
y = expression(paste(bold("Average # responded genes"))))
plot_PC2
# just selection
numResponded <- table(selResults$topo, selResults$responseToSel)
numResponded##
## FALSE TRUE
## ER 4647 32616
## BA 3494 35795
## WS 6519 30203numTopos <- as.vector(table(selResults$topo))
numResponded[, 2]/numTopos## ER BA WS
## 0.8752918 0.9110693 0.8224770#plot_numRespondedGenes <- ggplot(respondedGenes, aes(y = , x = topo)) +
# geom_violin(fill = genotypeColor, color = "black", trim = TRUE) +
# geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
# stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
# hjust = 1.25, vjust = 3, color = "black", aes(label = round(..y.., digits = 3))) +
# stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
# theme_pubclean() +
# theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
# axis.title.x = element_text(size=16, face="bold"),
# axis.title.y = element_text(size=16, face="bold"),
# axis.text.x = element_text(size=10, face="bold"),
# axis.text.y = element_text(size=10, face="bold")) +
# scale_fill_manual(values = c("sel" = genotypeColor, "neutr" = genotypeColor)) +
# labs(x = "Network topology",
# y = expression(paste(bold("Selective pressure "), "|", bold(p), "|")))
#plot_numRespondedGenes
#ggsave(filename = sprintf("plot_numRespondedGenes.png"),
# plot = plot_numRespondedGenes,
# device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
# dpi = 300, limitsize = TRUE)6.2 Average expression variance, gen 1
numPermutations = 10000
binNum = 30
# observed MI
obs_MI_PC1 <- calcInformation(netSelResults$ave_varP_1,
netSelResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netSelResults$ave_varP_1,
netSelResults$PC2, binNum)
# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)
for(permNum in 1:numPermutations)
{
MI_nullDistr_PC1[permNum] <-
calcInformation(netSelResults$ave_varP_1,
sample(netSelResults$PC1,
size = length(netSelResults$PC1)),
binNum)
MI_nullDistr_PC2[permNum] <-
calcInformation(netSelResults$ave_varP_1,
sample(netSelResults$PC2,
size = length(netSelResults$PC2)),
binNum)
}
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# title
if(pval_MI_absInStrT == 1e-04) {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ", round(pval_MI_absInStrT, digits = 2)))
}
if(pval_MI_absOutStrT == 1e-04) {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ", round(pval_MI_absOutStrT, digits = 2)))
}
# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_ave_varP_1-PCs.txt"))
cat(paste0("Variables: ave_varP_1; PC1\n",
"Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n",
"Variables: ave_varP_1; PC2\n",
"Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))## Variables: ave_varP_1; PC1
## Observed MI: 0.206465662538958; pval: 1e-04
## Variables: ave_varP_1; PC2
## Observed MI: 0.206738247502956; pval: 1e-04sink()lmModel <- lm(ave_varP_1 ~ PC1 + PC2,
data = netSelResults)
summary(lmModel)##
## Call:
## lm(formula = ave_varP_1 ~ PC1 + PC2, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -296.37 -70.77 -9.85 57.67 512.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 363.9316 1.8418 197.593 <2e-16 ***
## PC1 -6.1872 0.7399 -8.362 <2e-16 ***
## PC2 13.2576 1.0290 12.884 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 100.9 on 2997 degrees of freedom
## Multiple R-squared: 0.07298, Adjusted R-squared: 0.07236
## F-statistic: 118 on 2 and 2997 DF, p-value: < 2.2e-16r.squaredGLMM(lmModel)## R2m R2c
## [1,] 0.07293031 0.07293031# partial R^2
library(rstatix)
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]
# R2m for plotting
partial_R2m_absInStr_num <- round(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- round(partial_R2_PC2, digits = 2)
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))
# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 2)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 3)
# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)
# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
{pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else
{pval_coef2_title = paste0("p = ", pval_explVar2)}
coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))
plot_PC1 <- ggplot(netSelResults, aes(y = ave_varP_1, x = PC1)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
geom_abline(slope = lmModel$coefficients[2],
intercept = lmModel$coefficients[1],
color = "black", linetype = "dashed", size = 1) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_colour_manual(values = topoColors) +
labs(x = expression(bold("PC1 (diameter + centralization)")),
y = expression(paste(bold("Expression variance, gen. 1"))),
title = coef_pval_explVar1_title,
subtitle = R2m_absInStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -3, y = 750, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4, hjust = 0)
plot_PC1## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
plot_PC2 <- ggplot(netSelResults, aes(y = ave_varP_1, x = PC2)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
geom_abline(slope = lmModel$coefficients[3],
intercept = lmModel$coefficients[1],
color = "black", linetype = "dashed", size = 1) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_colour_manual(values = topoColors) +
labs(x = expression(bold("PC2 (average degree)")),
y = expression(paste(bold("Expression variance, gen. 1"))),
title = coef_pval_explVar2_title,
subtitle = R2m_absOutStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -3, y = 750, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4, hjust = 0)
plot_PC2## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
plot_netPropertiesFigure <- plot_grid(plot_PC1, plot_PC2,
labels = "AUTO",
label_size = 20,
label_fontface = "bold")## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'ggsave(filename = sprintf("plot_netPropertiesFigure_averageExprVar.png"),
plot = plot_netPropertiesFigure,
bg = "white",
path = pathToPlotsFolder,
device = "png",
scale = 2, height = 800, width = 1500, units = "px",
dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_netPropertiesFigure_averageExprVar.tiff"),
plot = plot_netPropertiesFigure,
bg = "white",
path = pathToPlotsFolder,
device = "tiff",
scale = 2, height = 800, width = 1500, units = "px",
dpi = 300, limitsize = TRUE)6.3 Change of expression variance
numPermutations = 10000
binNum = 30
# observed MI
obs_MI_PC1 <- calcInformation(netSelResults$ave_relDeltaVar_10000,
netSelResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netSelResults$ave_relDeltaVar_10000,
netSelResults$PC2, binNum)
# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)
for(permNum in 1:numPermutations)
{
MI_nullDistr_PC1[permNum] <-
calcInformation(netSelResults$ave_relDeltaVar_10000,
sample(netSelResults$PC1,
size = length(netSelResults$PC1)),
binNum)
MI_nullDistr_PC2[permNum] <-
calcInformation(netSelResults$ave_relDeltaVar_10000,
sample(netSelResults$PC2,
size = length(netSelResults$PC2)),
binNum)
}
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# title
if(pval_MI_absInStrT == 1e-04) {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ", round(pval_MI_absInStrT, digits = 2)))
}
if(pval_MI_absOutStrT == 1e-04) {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ", round(pval_MI_absOutStrT, digits = 2)))
}
# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_ave_relDeltaVar_10000-PCs.txt"))
cat(paste0("Variables: ave_relDeltaVar_10000; PC1\n",
"Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n",
"Variables: ave_relDeltaVar_10000; PC2\n",
"Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))## Variables: ave_relDeltaVar_10000; PC1
## Observed MI: 0.210243558120783; pval: 1e-04
## Variables: ave_relDeltaVar_10000; PC2
## Observed MI: 0.194051558163579; pval: 1e-04sink()lmModel <- lm(ave_relDeltaVar_10000 ~ PC1 + PC2,
data = netSelResults)
summary(lmModel)##
## Call:
## lm(formula = ave_relDeltaVar_10000 ~ PC1 + PC2, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04093 -0.01302 -0.00595 0.00257 0.60148
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9184551 0.0006813 -1348.061 < 2e-16 ***
## PC1 0.0023709 0.0002719 8.718 < 2e-16 ***
## PC2 -0.0010134 0.0003826 -2.649 0.00813 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0367 on 2900 degrees of freedom
## (97 observations deleted due to missingness)
## Multiple R-squared: 0.02782, Adjusted R-squared: 0.02715
## F-statistic: 41.5 on 2 and 2900 DF, p-value: < 2.2e-16r.squaredGLMM(lmModel)## R2m R2c
## [1,] 0.0278045 0.0278045# partial R^2
library(rstatix)
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]
# R2m for plotting
partial_R2m_absInStr_num <- round(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- round(partial_R2_PC2, digits = 2)
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))
# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 2)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 3)
#if(coef_explVar2 == -0.0096){coef_explVar2 = -0.009}
# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)
# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
{pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else
{pval_coef2_title = paste0("p = ", pval_explVar2)}
coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))
plot_PC1 <- ggplot(netSelResults, aes(y = ave_relDeltaVar_10000, x = PC1)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
geom_abline(slope = lmModel$coefficients[2],
intercept = lmModel$coefficients[1],
color = "black", linetype = "dashed", size = 1) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_colour_manual(values = topoColors) +
labs(x = expression(bold("PC1 (diameter + centralization)")),
y = expression(bold("Rel."~Delta~"expr. variance")),
title = coef_pval_explVar1_title,
subtitle = R2m_absInStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -4, y = -0.4, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4, hjust = 0)
plot_PC1## Warning: Removed 97 rows containing missing values (geom_point).## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
plot_PC2 <- ggplot(netSelResults, aes(y = ave_relDeltaVar_10000, x = PC2)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
geom_abline(slope = lmModel$coefficients[3],
intercept = lmModel$coefficients[1],
color = "black", linetype = "dashed", size = 1) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "right") +
scale_colour_manual(values = topoColors) +
labs(x = expression(bold("PC2 (average degree)")),
y = expression(bold("Rel."~Delta~"expr. variance")),
title = coef_pval_explVar2_title,
subtitle = R2m_absOutStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -3, y = -0.4, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4, hjust = 0)
plot_PC2## Warning: Removed 97 rows containing missing values (geom_point).
## Warning: is.na() applied to non-(list or vector) of type 'expression'
plot_netPropertiesFigure <- plot_grid(plot_PC1, plot_PC2,
labels = "AUTO",
label_size = 20,
label_fontface = "bold")## Warning: Removed 97 rows containing missing values (geom_point).
## Warning: is.na() applied to non-(list or vector) of type 'expression'## Warning: Removed 97 rows containing missing values (geom_point).## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'ggsave(filename = sprintf("plot_netPropertiesFigure_averageRelDeltaVar.png"),
plot = plot_netPropertiesFigure,
bg = "white",
path = pathToPlotsFolder,
device = "png",
scale = 2, height = 800, width = 1500, units = "px",
dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_netPropertiesFigure_averageRelDeltaVar.tiff"),
plot = plot_netPropertiesFigure,
bg = "white",
path = pathToPlotsFolder,
device = "tiff",
scale = 2, height = 800, width = 1500, units = "px",
dpi = 300, limitsize = TRUE)lmeModel <- lm(ave_s_g_area_abs ~ PC1 + PC2,
data = netSelResults)
summary(lmeModel)##
## Call:
## lm(formula = ave_s_g_area_abs ~ PC1 + PC2, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64402 -0.02321 0.01087 0.03897 0.12219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7741980 0.0011872 652.125 < 2e-16 ***
## PC1 -0.0031471 0.0004769 -6.599 4.88e-11 ***
## PC2 -0.0095611 0.0006633 -14.415 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06503 on 2997 degrees of freedom
## Multiple R-squared: 0.07737, Adjusted R-squared: 0.07676
## F-statistic: 125.7 on 2 and 2997 DF, p-value: < 2.2e-16r.squaredGLMM(lmeModel)## R2m R2c
## [1,] 0.07732693 0.077326936.4 Selective pressure
6.4.1 Neutrality
numPermutations = 10000
binNum = 30
# observed MI
obs_MI_PC1 <- calcInformation(netNeuResults$ave_s_g_area_abs,
netNeuResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netNeuResults$ave_s_g_area_abs,
netNeuResults$PC2, binNum)
# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)
for(permNum in 1:numPermutations)
{
MI_nullDistr_PC1[permNum] <-
calcInformation(netNeuResults$ave_s_g_area_abs,
sample(netNeuResults$PC1,
size = length(netNeuResults$PC1)),
binNum)
MI_nullDistr_PC2[permNum] <-
calcInformation(netNeuResults$ave_s_g_area_abs,
sample(netNeuResults$PC2,
size = length(netNeuResults$PC2)),
binNum)
}
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# title
if(pval_MI_absInStrT == 1e-04) {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ", round(pval_MI_absInStrT, digits = 2)))
}
if(pval_MI_absOutStrT == 1e-04) {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ", round(pval_MI_absOutStrT, digits = 2)))
}
# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_s_g_area_abs-PCs_neutrality.txt"))
cat(paste0("Variables: ave_s_g_area_abs; PC1\n",
"Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n",
"Variables: ave_s_g_area_abs; PC2\n",
"Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))## Variables: ave_s_g_area_abs; PC1
## Observed MI: 0.147939366727152; pval: 0.7226
## Variables: ave_s_g_area_abs; PC2
## Observed MI: 0.150486308414081; pval: 0.5815sink()
plot_hist_MI_obs_PC1 <- ggplot(data = data.frame(MI = MI_nullDistr_PC1), aes(x = MI)) +
geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
geom_vline(xintercept = obs_MI_PC1, col = MIColor, lwd = 2, lty = 2) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
labs(x = expression(paste(bold("Mutual Information (PC1)"))), y = "count",
title = MI_obs_explVar1_title_with_pval) +
annotate('text', x = obs_MI_PC1, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC1
ggsave(filename = sprintf("plot_hist_MI_obs_PC1_neu.png"),
path = pathToPlotsFolder,
plot = plot_hist_MI_obs_PC1,
device = "png", scale = 2, width = 6, height = 6, units = "cm",
dpi = 300, limitsize = TRUE)
plot_hist_MI_obs_PC2 <- ggplot(data = data.frame(MI = MI_nullDistr_PC2), aes(x = MI)) +
geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
geom_vline(xintercept = obs_MI_PC2, col = MIColor, lwd = 2, lty = 2) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
labs(x = expression(paste(bold("Mutual Information (PC2)"))), y = "count",
title = MI_obs_explVar2_title_with_pval) +
annotate('text', x = obs_MI_PC2, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC2
ggsave(filename = sprintf("plot_hist_MI_obs_PC2_neu.png"),
path = pathToPlotsFolder,
plot = plot_hist_MI_obs_PC2,
device = "png", scale = 2, width = 6, height = 6, units = "cm",
dpi = 300, limitsize = TRUE)indeptest <- function(model) {
return(Box.test(resid(model)[order(fitted(model))], type = "Ljung-Box"))
}
lmModel <- lm(ave_s_g_area_abs ~ PC1 + PC2,
data = netNeuResults)
summary(lmModel)##
## Call:
## lm(formula = ave_s_g_area_abs ~ PC1 + PC2, data = netNeuResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.913e-03 -3.902e-04 -1.172e-05 3.846e-04 1.936e-03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.078e-03 1.065e-05 664.782 <2e-16 ***
## PC1 -2.836e-07 4.277e-06 -0.066 0.947
## PC2 -1.022e-07 5.948e-06 -0.017 0.986
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0005832 on 2997 degrees of freedom
## Multiple R-squared: 1.566e-06, Adjusted R-squared: -0.0006658
## F-statistic: 0.002346 on 2 and 2997 DF, p-value: 0.9977r.squaredGLMM(lmModel)## R2m R2c
## [1,] 1.564803e-06 1.564803e-06vif(lmModel)## PC1 PC2
## 1 1shapiro.test(lmModel[['residuals']])##
## Shapiro-Wilk normality test
##
## data: lmModel[["residuals"]]
## W = 0.99936, p-value = 0.4045indeptest(lmModel)##
## Box-Ljung test
##
## data: resid(model)[order(fitted(model))]
## X-squared = 0.81939, df = 1, p-value = 0.3654partial_eta_squared(lmModel)## PC1 PC2
## 1.467371e-06 9.847603e-08# partial R^2
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]
# R2m for plotting
partial_R2m_absInStr_num <- signif(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- signif(partial_R2_PC2, digits = 2)
if(partial_R2m_absInStr_num == 1.5e-06){partial_R2m_absInStr_num = "1.5 x 10^{-6}"}
if(partial_R2m_absOutStr_num == 9.8e-08){partial_R2m_absOutStr_num = "9.8 x 10^{-8}"}
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))
# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 2)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 2)
if(coef_explVar1 == -2.8e-07){coef_explVar1 = "-2.8 x 10^{-7}"}
if(coef_explVar2 == -1e-07){coef_explVar2 = "-1 x 10^{-7}"}
# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)
# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
{pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else
{pval_coef2_title = paste0("p = ", pval_explVar2)}
coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))
plot_constVar_resid <- ggplot(data = data.frame("Fitted_values" = fitted(lmModel),
"Pearsons_residuals" = resid(lmModel, type = "pearson")),
aes(x = Fitted_values, y = Pearsons_residuals)) +
geom_point(alpha = 0.2, size = 0.5) +
theme_bw() +
theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
axis.title.x = element_text(size=8, face="bold"),
axis.title.y = element_text(size=8, face="bold"),
axis.text.x = element_text(size=6, face="bold"),
axis.text.y = element_text(size=6, face="bold")) +
annotate("label", x = 0.0070775, y = 0.0015, label = "Neutrality", size = 3) +
labs(x = "Fitted values", y = "Pearson's residuals")
plot_constVar_qqResid <- ggplot(data = data.frame("Pearsons_residuals" = resid(lmModel, type = "pearson")),
aes(sample = Pearsons_residuals)) +
stat_qq(alpha = 0.2, size = 0.5) + stat_qq_line() +
theme_bw() +
theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
axis.title.x = element_text(size=8, face="bold"),
axis.title.y = element_text(size=8, face="bold"),
axis.text.x = element_text(size=6, face="bold"),
axis.text.y = element_text(size=6, face="bold")) +
labs(x = "Theoretical quantiles", y = "Sample quantiles",
title = "Normal Q-Q plot, residuals")
plot_constVar_resid
plot_constVar_qqResid
plot_PC1_selpress_neu <- ggplot(netNeuResults, aes(y = ave_s_g_area_abs, x = PC1)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "bottom",
legend.title = element_text(size=12),
legend.text = element_text(size=12)) +
scale_colour_manual(values = topoColors,
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15),
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
labs(x = expression(bold("PC1 (diameter + centralization)")),
y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
title = coef_pval_explVar1_title,
subtitle = R2m_absInStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -4, y = 0.95, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4, hjust = 0) +
ylim(0, 1)
plot_PC1_selpress_neu## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
plot_PC2_selpress_neu <- ggplot(netNeuResults, aes(y = ave_s_g_area_abs, x = PC2)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "bottom",
legend.title = element_text(size=12),
legend.text = element_text(size=12)) +
scale_colour_manual(values = topoColors,
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15),
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
labs(x = expression(bold("PC2 (average degree)")),
y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
title = coef_pval_explVar2_title,
subtitle = R2m_absOutStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -3, y = 0.95, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4, hjust = 0) +
scale_x_continuous(n.breaks = 4) +
ylim(0, 1)
plot_PC2_selpress_neu## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
6.4.2 Selection
numPermutations = 10000
binNum = 30
# observed MI
obs_MI_PC1 <- calcInformation(netSelResults$ave_s_g_area_abs,
netSelResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netSelResults$ave_s_g_area_abs,
netSelResults$PC2, binNum)
# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)
for(permNum in 1:numPermutations)
{
MI_nullDistr_PC1[permNum] <-
calcInformation(netSelResults$ave_s_g_area_abs,
sample(netSelResults$PC1,
size = length(netSelResults$PC1)),
binNum)
MI_nullDistr_PC2[permNum] <-
calcInformation(netSelResults$ave_s_g_area_abs,
sample(netSelResults$PC2,
size = length(netSelResults$PC2)),
binNum)
}
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# pvals for observed MI and R
pval_MI_absInStrT <-
(length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
(length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations
# title
if(pval_MI_absInStrT == 1e-04) {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar1_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ", round(pval_MI_absInStrT, digits = 2)))
}
if(pval_MI_absOutStrT == 1e-04) {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
MI_obs_explVar2_title_with_pval <-
TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ", round(pval_MI_absOutStrT, digits = 2)))
}
# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_s_g_area_abs-PCs.txt"))
cat(paste0("Variables: ave_s_g_area_abs; PC1\n",
"Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n",
"Variables: ave_s_g_area_abs; PC2\n",
"Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))## Variables: ave_s_g_area_abs; PC1
## Observed MI: 0.265906498613406; pval: 1e-04
## Variables: ave_s_g_area_abs; PC2
## Observed MI: 0.263734712146046; pval: 1e-04sink()
plot_hist_MI_obs_PC1 <- ggplot(data = data.frame(MI = MI_nullDistr_PC1), aes(x = MI)) +
geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
geom_vline(xintercept = obs_MI_PC1, col = MIColor, lwd = 2, lty = 2) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
labs(x = expression(paste(bold("Mutual Information (PC1)"))), y = "count",
title = MI_obs_explVar1_title_with_pval) +
annotate('text', x = obs_MI_PC1, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC1
ggsave(filename = sprintf("plot_hist_MI_obs_PC1_sel.png"),
path = pathToPlotsFolder,
plot = plot_hist_MI_obs_PC1,
device = "png", scale = 2, width = 6, height = 6, units = "cm",
dpi = 300, limitsize = TRUE)
plot_hist_MI_obs_PC2 <- ggplot(data = data.frame(MI = MI_nullDistr_PC2), aes(x = MI)) +
geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
geom_vline(xintercept = obs_MI_PC2, col = MIColor, lwd = 2, lty = 2) +
theme_pubclean() +
theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold")) +
labs(x = expression(paste(bold("Mutual Information (PC2)"))), y = "count",
title = MI_obs_explVar2_title_with_pval) +
annotate('text', x = obs_MI_PC2, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC2
ggsave(filename = sprintf("plot_hist_MI_obs_PC2_sel.png"),
path = pathToPlotsFolder,
plot = plot_hist_MI_obs_PC2,
device = "png", scale = 2, width = 6, height = 6, units = "cm",
dpi = 300, limitsize = TRUE)indeptest <- function(model) {
return(Box.test(resid(model)[order(fitted(model))], type = "Ljung-Box"))
}
lmModel <- lm(ave_s_g_area_abs ~ PC1 + PC2,
data = netSelResults)
summary(lmModel)##
## Call:
## lm(formula = ave_s_g_area_abs ~ PC1 + PC2, data = netSelResults)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64402 -0.02321 0.01087 0.03897 0.12219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7741980 0.0011872 652.125 < 2e-16 ***
## PC1 -0.0031471 0.0004769 -6.599 4.88e-11 ***
## PC2 -0.0095611 0.0006633 -14.415 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06503 on 2997 degrees of freedom
## Multiple R-squared: 0.07737, Adjusted R-squared: 0.07676
## F-statistic: 125.7 on 2 and 2997 DF, p-value: < 2.2e-16r.squaredGLMM(lmModel)## R2m R2c
## [1,] 0.07732693 0.07732693vif(lmModel)## PC1 PC2
## 1 1shapiro.test(lmModel[['residuals']])##
## Shapiro-Wilk normality test
##
## data: lmModel[["residuals"]]
## W = 0.81175, p-value < 2.2e-16indeptest(lmModel)##
## Box-Ljung test
##
## data: resid(model)[order(fitted(model))]
## X-squared = 0.91059, df = 1, p-value = 0.34partial_eta_squared(lmModel)## PC1 PC2
## 0.01432150 0.06483836# partial R^2
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]
# R2m for plotting
partial_R2m_absInStr_num <- round(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- round(partial_R2_PC2, digits = 2)
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))
# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 1)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 2)
if(coef_explVar2 == -0.0096){coef_explVar2 = -0.009}
# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)
# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
{pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else
{pval_coef2_title = paste0("p = ", pval_explVar2)}
if(pval_explVar1 == 4.9e-11) {pval_coef1_title = paste0("p = 4.9 x 10^{-11}")}
coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))
plot_constVar_resid_sel <- ggplot(data = data.frame("Fitted_values" = fitted(lmModel),
"Pearsons_residuals" = resid(lmModel, type = "pearson")),
aes(x = Fitted_values, y = Pearsons_residuals)) +
geom_point(alpha = 0.2, size = 0.5) +
theme_bw() +
theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
axis.title.x = element_text(size=8, face="bold"),
axis.title.y = element_text(size=8, face="bold"),
axis.text.x = element_text(size=6, face="bold"),
axis.text.y = element_text(size=6, face="bold")) +
annotate("label", x = 0.75, y = -0.5, label = "Selection", size = 3) +
labs(x = "Fitted values", y = "Pearson's residuals")
plot_constVar_qqResid_sel <- ggplot(data = data.frame("Pearsons_residuals" = resid(lmModel, type = "pearson")),
aes(sample = Pearsons_residuals)) +
stat_qq(alpha = 0.2, size = 0.5) + stat_qq_line() +
theme_bw() +
theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
axis.title.x = element_text(size=8, face="bold"),
axis.title.y = element_text(size=8, face="bold"),
axis.text.x = element_text(size=6, face="bold"),
axis.text.y = element_text(size=6, face="bold")) +
labs(x = "Theoretical quantiles", y = "Sample quantiles",
title = "Normal Q-Q plot, residuals")
plot_constVar_resid_sel
plot_constVar_qqResid_sel
plot_ModelDiagnostics <- plot_grid(plot_constVar_resid_sel, plot_constVar_qqResid_sel,
plot_constVar_resid, plot_constVar_qqResid,
ncol = 2,
labels = "AUTO")
# save to plots folder
ggsave(filename = sprintf("plot_ModelDiagnostics_networkProperties.png"),
path = pathToPlotsFolder,
plot = plot_ModelDiagnostics,
device = "png", scale = 1.2, width = 12, height = 12, units = "cm",
dpi = 300, limitsize = TRUE,
bg = "white")plot_PC1_selpress_sel <- ggplot(netSelResults, aes(y = ave_s_g_area_abs, x = PC1)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2) +
geom_quantile(quantiles = c(.5), color = "black", size = 0.75) +
geom_quantile(quantiles = c(.25, .75), color = "black", size = 0.5, linetype = 2) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "bottom",
legend.title = element_text(size=12),
legend.text = element_text(size=12)) +
scale_colour_manual(values = topoColors,
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15),
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
labs(x = expression(bold("PC1 (diameter + centralization)")),
y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
title = coef_pval_explVar1_title,
subtitle = R2m_absInStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -4, y = 0.95, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4, hjust = 0) +
ylim(0, 1)
plot_PC1_selpress_sel## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
plot_PC2_selpress_sel <- ggplot(netSelResults, aes(y = ave_s_g_area_abs, x = PC2)) +
geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2) +
geom_quantile(quantiles = c(.5), color = "black", size = 0.75) +
geom_quantile(quantiles = c(.25, .75), color = "black", size = 0.5, linetype = 2) +
theme_bw() +
theme(plot.title = element_text(size=12, face="bold",),
plot.subtitle = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(size=10, face="bold"),
axis.text.y = element_text(size=10, face="bold"),
legend.position = "bottom",
legend.title = element_text(size=12),
legend.text = element_text(size=12)) +
scale_colour_manual(values = topoColors,
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15),
labels = c("random (Erdős–Rényi)",
"scale-free (Barabási–Albert)",
"small-world (Watts–Strogatz)")) +
labs(x = expression(bold("PC2 (average degree)")),
y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
title = coef_pval_explVar2_title,
subtitle = R2m_absOutStr_text,
color = "Topology",
shape = "Topology") +
annotate('text', x = -3, y = 0.95, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4, hjust = 0) +
scale_x_continuous(n.breaks = 4) +
ylim(0, 1)
plot_PC2_selpress_sel## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
7 Final plot
jointTitle_sel <- ggdraw() + draw_label("Selection",
size = 20,
fontface = 'bold')
jointTitle_neu <- ggdraw() + draw_label("Neutrality",
size = 20,
fontface = 'bold')
jointTitle_combined <- cowplot::plot_grid(NULL, jointTitle_sel, NULL,
NULL, jointTitle_neu, NULL,
labels = c("", "", "", "", "", ""),
ncol = 6,
rel_widths = c(0.5, 1, 0.5, 0.5, 1, 0.5))
plot_netPropertiesFigure_body <- ggpubr::ggarrange(plot_PC1_selpress_sel, plot_PC2_selpress_sel,
plot_PC1_selpress_neu, plot_PC2_selpress_neu,
labels = "AUTO", font.label = list(size = 20, face = "bold"),
ncol = 4, nrow = 1,
common.legend = TRUE, legend = "bottom")## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'plot_netPropertiesFigure <- cowplot::plot_grid(jointTitle_combined,
plot_netPropertiesFigure_body,
ncol = 1,
rel_heights = c(0.1, 1))
ggsave(filename = sprintf("plot_netPropertiesFigure.png"),
plot = plot_netPropertiesFigure,
bg = "white",
path = pathToPlotsFolder,
device = "png",
scale = 2.1, height = 800, width = 2250, units = "px",
dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_netPropertiesFigure.tiff"),
plot = plot_netPropertiesFigure,
bg = "white",
path = pathToPlotsFolder,
device = "tiff",
scale = 2.1, height = 900, width = 2250, units = "px",
dpi = 300, limitsize = TRUE)sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] webshot_0.5.2 htmltools_0.5.2 formattable_0.2.1 dplyr_1.0.7
## [5] rstatix_0.7.0 FSA_0.9.3 factoextra_1.0.7 ade4_1.7-15
## [9] corrplot_0.90 Hmisc_4.3-1 Formula_1.2-3 survival_3.2-7
## [13] lattice_0.20-38 reshape2_1.4.3 latex2exp_0.4.0 RColorBrewer_1.1-2
## [17] car_3.0-11 carData_3.0-4 lme4_1.1-27.1 Matrix_1.2-18
## [21] infotheo_1.2.0 cowplot_1.1.1 gridExtra_2.3 ggridges_0.5.2
## [25] ggpubr_0.4.0 ggplot2_3.3.5 MuMIn_1.43.17 nlme_3.1-144
## [29] rmarkdown_2.10
##
## loaded via a namespace (and not attached):
## [1] matrixStats_0.60.0 tools_3.6.3 backports_1.2.1
## [4] bslib_0.2.5.1 utf8_1.2.2 R6_2.5.1
## [7] rpart_4.1-15 DBI_1.1.0 colorspace_2.0-2
## [10] nnet_7.3-12 withr_2.4.2 tidyselect_1.1.1
## [13] curl_4.3.2 compiler_3.6.3 quantreg_5.86
## [16] htmlTable_1.13.3 SparseM_1.81 labeling_0.4.2
## [19] sass_0.4.0 scales_1.1.1 checkmate_2.0.0
## [22] stringr_1.4.0 digest_0.6.28 foreign_0.8-75
## [25] minqa_1.2.4 rio_0.5.27 base64enc_0.1-3
## [28] jpeg_0.1-8.1 pkgconfig_2.0.3 dunn.test_1.3.5
## [31] highr_0.9 fastmap_1.1.0 htmlwidgets_1.5.3
## [34] rlang_0.4.12 readxl_1.3.1 rstudioapi_0.13
## [37] farver_2.1.0 jquerylib_0.1.4 generics_0.1.0
## [40] jsonlite_1.7.2 acepack_1.4.1 zip_2.2.0
## [43] magrittr_2.0.1 Rcpp_1.0.7 munsell_0.5.0
## [46] fansi_0.5.0 abind_1.4-5 lifecycle_1.0.1
## [49] stringi_1.7.3 yaml_2.2.1 MASS_7.3-57
## [52] plyr_1.8.6 grid_3.6.3 ggrepel_0.9.1
## [55] forcats_0.5.1 crayon_1.4.2 haven_2.4.3
## [58] splines_3.6.3 hms_1.1.0 knitr_1.33
## [61] pillar_1.6.4 boot_1.3-25 ggsignif_0.6.2
## [64] stats4_3.6.3 glue_1.5.0 evaluate_0.14
## [67] latticeExtra_0.6-29 data.table_1.14.0 vctrs_0.3.8
## [70] png_0.1-7 nloptr_1.2.2.2 MatrixModels_0.5-0
## [73] cellranger_1.1.0 gtable_0.3.0 purrr_0.3.4
## [76] tidyr_1.1.3 xfun_0.25 openxlsx_4.2.4
## [79] broom_0.7.9 conquer_1.0.2 tibble_3.1.6
## [82] cluster_2.1.0 ellipsis_0.3.2packageVersion('igraph')## [1] '1.2.4.2'packageVersion('intergraph')## [1] '2.0.2'packageVersion('lme4')## [1] '1.1.27.1'packageVersion('nlme')## [1] '3.1.144'packageVersion('MuMIn')## [1] '1.43.17'packageVersion('infotheo')## [1] '1.2.0'packageVersion('car')## [1] '3.0.11'packageVersion('ade4')## [1] '1.7.15'